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Abstract 


We report on integral field unit (IFU) measurements of the host of the radio source 4C+37.11. This massive 
elliptical contains the only resolved double compact nucleus at parsec-scale separation, likely a bound 
supermassive black hole binary (SMBHB). i-band photometry and GMOS-N IFU spectroscopy show that the 
galaxy has a large rj = 1"5 core and that the stellar velocity dispersion increases inside of a radius of influence 


soi © 1"3. Jeans Anisotropic Modeling analysis of the core infers a total SMBHB mass of 2.8103 


x 10M, 


making this one of the most massive black hole systems known. Our data indicate that there has been significant 


scouring of the central kiloparsec of the host galaxy. 


Unified Astronomy Thesaurus concepts: Supermassive black holes (1663); Giant elliptical galaxies (651); Brightest 


cluster galaxies (181) 


Supporting material: data behind figure 


1. Introduction 


Maness et al. (2004) found that the radio galaxy 4C+37.11 
(04024-379) contains two central, compact, flat-spectrum, 
variable components (designated C1 and C2) with a Very 
Long Baseline Array (VLBA)-measured separation 7.3 pc at 
z=0.055 and argued that this is a gravitationally bound 
supermassive black hole binary (SMBHB). Subsequent VLBA 
observations by Rodriguez et al. (2006) have bolstered that 
claim and even detected a possible relative proper motion of the 
binary components (Bansal et al. 2017). The binary's host is 
remarkable, being an extremely massive elliptical embedded in 
a bright X-ray halo (Romani et al. 2014; Andrade-Santos et al. 
2016), which represents a cluster's worth of stars and mass in a 
relatively isolated galaxy. Such “fossil clusters" are believed to 
be the product of multiple major mergers. We report here on a 
kinematic study of the central region of this galaxy, which 
supports a scenario where the source contains a supermassive 
black hole pair whose present (and former) interactions have 
scoured the galaxy core. The estimated SMBHB mass is very 
large, the second largest kinematically measured value in the 
local Universe. 

We assume a 2H9—69.6, Qy=0.286, and 04,—0.714 
cosmology (Bennett et al. 2014), so at the redshift of z — 0.055, 
4C--37.11 has a luminosity distance D, = 246.9 Mpc, angular 
diameter distance D4-— 221.9 Mpc, and a projected scale of 
1.076 kpc! 


2. Observations, Data Reduction, and Measurement 


Our kinematic study uses two observations. A 2008 
December 12 WIYN 3.6m 2200s OPTIC i band image 
presented in Romani et al. (2014) has image quality FWHM 
<0”6 and measures the galaxy surface brightness profile. For 


Original content from this work may be used under the terms 

BY of the Creative Commons Attribution 4.0 licence. Any further 
distribution of this work must maintain attribution to the author(s) and the title 
of the work, journal citation and DOI. 


spectroscopy, we use a 2015 December 1 Gemini 8.1 m 
6 x 1200s GMOS-N integral field unit (IFU) observation of 
the galaxy center, downloaded from the Gemini Science 
archive. These observations used the B600 grating, covering 
6300-9195 Å at dispersion 0.47 A pixel" !. Conditions were 
good, with low airmass throughout. The IFU feeds a fiber array 
with a 0"22 pitch, and with the observations dithered on the 
sky, IRAF GMOS pipeline processing calibrated and cleaned 
the spectra, assembling a data cube covering 4"5 x 6"2 at 0” 1 
scale, centered on the active galactic nucleus (AGN), with 
reduced net exposure at the edges. Measurements of images of 
the point-spread function (PSF)-distributed line wings of the 
AGN core find an angular resolution of 0764 (FWHM) and 
sky lines give a median effective spectral resolution of 
2.0 AFWHM. 

For our kinematic study, we fit the emission line-free range 
7810-9190 A, masking a small range near 8230 A affected by a 
chip gap. We log rebin the spectra over this range to a velocity 
scale of 15.3kms . Reliable measurements of the line-of- 
sight velocity distribution (LOSVD) require a signal-to-noise 
ratio (S/N) of —20 per resolution element. Our data cube 
provides only S/N ~7 per spaxel in the galaxy center, falling 
to S/N ~1 in the more poorly exposed outskirts of the data 
cube. Thus, we use Voronoi binning (Cappellari & Copin 2003) 
to combine the spaxels. With a target S/N ~20 at the rebinned 
velocity scale, we get 31 bins (Figure 1); edge bins generally 
do not reach the S/N target. Imperfectly excised cosmic-ray 
events and poor sky subtraction affect many of the edge pixels 
with few IFU exposures. We therefore used iterative sigma 
clipping to excise outlier pixels from each wavelength in the 
Voronoi tiles. The means of the clipped tile wavelength 
measurements form the final spectra, with their dispersion as 
the associated error vector. Velocity structure is measured with 
the Penalized PiXel Fitting method (pPXF) of Cappellari 
(2023), which convolves template stellar spectra with a 
Gauss-Hermite decomposition of the LOSVD. With the 
limited S/N, we fit only V and c, skipping the higher order 
h, and h4 moments. 
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Figure 1. Radial velocity (V) and velocity dispersion (c) maps of pPXF measurements in the Voronoi bins. The origin is the AGN location, determined from the wings of the 


[N N] emission lines. Solid contours mark 0.9x,0.8x,.. 


For the stellar templates, we used a combination of GO-M2 
giants from the Indo-US Stellar Library (Valdes et al. 2004). 
These templates have a slightly finer scale of O. 44 A pixel !, 
with a median spectral resolution of ~1 A. To match the 
templates to the host galaxy spectra, we convolve the templates 
using the difference in the spectral resolutions, following 
Section 2.2 in Cappellari (2017). pPXF uses additive and 
multiplicative polynomials in the fitting model for both the 
target and template spectra to reduce sensitivity to imperfect 
continuum fluxing and veiling components. We varied the 
order of these polynomials, finding that both additive and 
multiplicative terms required order >4 to achieve stable fits. 
Here we use fourth-order additive and sixth-order multi- 
plicative Legendre polynomials (4A+6M) as fiducial. Sample 
spectral fits at three different radii are shown in Figure 2. 

Radial velocity and velocity dispersion maps from the 
Voronoi tiles are shown in Figure 1. We detect little if any 
coherent rotation, typical for massive ellipticals. However, 
given the large size of the Voronoi bins, we cannot exclude a 
kinematically distinct core, sometimes seen among nonregular 
rotators (Cappellari 2016). We do see a substantial Keplerian- 
like peak in o(r), which flattens out at ~300 km s~! beyond 
~2". Although we lack kinematic data at the half-light radius 
(Re~ 12"; see below in the discussion of the Multi-Gaussian 
Expansion), we can follow Saglia et al. (2016, their Appendix 
A) and estimate the effective velocity dispersion by integrating 
the surface-brightness-weighted V,,, to our outermost bin, 
finding c, = 332+ 11 kms -! This likely overestimates e, as 
over a third of our bins are affected by the central velocity cusp. 

For dynamical modeling, we also require the host surface 
brightness profile. This was measured from the WIYN i’ image 
using the Multi-Gaussian Expansion (MGE) method of 
Cappellari (2002). We first masked stars and galaxies near 
4C+37.11. We then obtained the image PSF from 15 
unsaturated stars near the host, using the photutils effective 
PSF (ePSF) routine (Bradley et al. 2023). In host fitting we 
assume axisymmetry, with the surface brightness modeled as a 
sum of elliptical Gaussian components aligned along the 
direction of the photometric major axis (x^): 


Xx! WES Dod. 12 y? 1 
y= Y Sexp| -x+ (1) 


= Oj q; 


Here we have N Gaussian components with surface brightness S; 
and axial ratio qj. The position angles of the Gaussians were fixed 


.0.2x the peak surface brightness and dashed lines separate quadrants along the major and minor axes. 


0.0121 / 70.16" 


0.010 4 


0.008 4 


Relative Flux 


0.006 4 


00051  .. - em "UE MS oe ee 
0.0000 ] Fa ae Ein seh AA re ae TE 
—0.0025 . à i ri 


T T T T T T T 


r=1.41" 


0.006 4 
0.005 4 
0.004 4 
vial 

0.0025 l l l ul 

0.0000 |  sefanaie Ann eut et Nat Set rei SR PN ai pee oq Mt INS 

-0.0025 1 


Relative Flux 


Relative Flux 


0.001 4 


0.0025 4 
0.0000 4 
-0.0025 | 


T T T 
8000 8200 8400 8600 


Rest Frame Wavelength (A) 


T T -— 
7400 7600 7800 


Figure 2. Sample spectral fits (red) and residuals for S/N ~20 Voronoi bins at 
three different radii. For plotting purposes, the spectra were rebinned to 
61.2 km s^!, 4x the width used for the kinematic extraction. Grey bars show 
masked pixels from chip gaps and from 2c iterative clipping. 


at PAwgg = 7274 and qj; had a flat prior 0.75—1. Before fitting we 
convolve Equation (1) with an MGE expansion of the WIYN PSF. 
The MGE fits for the PSF and the i’ band image are summarized in 
Tables 1 and 2, respectively. Figure 3 shows the MGE fits for the 
IFU region and the entire region used in mass modeling. While 
somewhat triaxial at large radii, the axisymmetric assumption is 
adequate over the IFU kinematics region. Using the MGE, with the 
routine nge half light isophote of the package JamPy 
of Cappellari (2008) we find a half-light radius for 4C+-37.11 to be 
R,— 12"1. 

It was noted in Romani et al. (2014) that the host-surface- 
brightness profile flattened within ~1.5 kpc from an outer n = 4 
Sérsic profile, indicating that 4C+37.11 is a cored elliptical. 
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Figure 3. i’ image isophote contours (black) and MGE fit ellipses (red). Top: 
the full region. Bottom: the IFU region with kinematic measurements. The 
yellow regions correspond to masked pixels in the MGE fits. 


Table 1 
WIYN PSF MGE 
J Sj gj (arcsec) dj 
1 0.23 0.15 0.71 
2 1.38 0.27 0.93 
3 0.362 0.36 1.00 
4 0.10 0.56 0.96 
3 0.02 0.94 1.00 
6 0.002 21 1.00 


We measure the core radius r, by fitting a core-Sérsic model, 
with a point-source excess for the AGN light, to the i’ band 
image, using the reconstructed PSF for model convolution and 
minimizing with the 2D fitting software imfit (Erwin 2015). 
The core-Sérsic profile is described as 


2 o Nds 
I(r) = ZE + (2) | exp [t , Q) 
r Fe 


where 


I! = 27^ exp [by QV?rj/r,) /"], (3) 


Surti et al. 


Table 2 
PSF-convolved 4C4-37.11 MGE 
j log(Sj) (Lo pc?) 0; (arcsec) dj 
1 4.18 0.11 0.75 
2 3.18 1.25 0.75 
3 2.83 3.02 0.76 
4 2.29 6.38 0.75 
5 2.01 17.8 0.75 
Table 3 
2D Core-Sérsic + Point-source Model-fit Parameters 
Core Sérsic Parameters Value 
PA (deg) 81.9+0.3 
€ 0.293 + 0.003 
Sy. (mag arcsec ?) 18.26 + 0.03 
n 3.92 + 0.08 
r, (arcsec) 19.0 + 0.5 
ry (arcsec) 1.49 + 0.05 
y 0.40 0.03 
Point-source Parameters Value 
itor (mag) 20.2 + 0.3 
Note. 
Magnitudes $5, ; and i’ are corrected for an estimated A; = 2.2 mag extinction. 


as per Graham et al. (2003). For our fit, we fix the sharpness 
parameter at a = 100 (Graham et al. 2003; Trujillo et al. 2004) 
and derive the fit parameters in Table 3, where we have 
converted intensity at the break radius to a surface brightness 
Sp,’ We do indeed see that the outer galaxy has n ~ 4 and find 
a core radius r, = 1.60 + 0.05 kpc. For o > 10, we find that r, 
varies by «2o from this value. 


3. Dynamical Modeling 


Given the low S/N kinematics data extending only to ~3”, 
we do not attempt full Schwarzschild orbit-superposition 
modeling. Instead, we conducted axisymmetric Jeans 
Anisotropic Modeling (JAM) using the software JamPy 
(Cappellari 2008), which does not require large-radii 
kinematics to achieve good constraints (Simon et al. 2024) 
and produces black hole masses consistent with Schwarzschild 
analyses (Cappellari et al. 2010). 

Unlike the Schwarzschild method, axisymmetric JAM makes a 
number of assumptions about the distribution of stellar velocities. 
In particular, the velocity ellipsoid is either assumed to be aligned 
with spherical (r, $, 0) or cylindrical (R, ġ, z) coordinates, neither 
of which can fully represent a real galaxy. Spherically aligned 
velocity ellipsoids are typically good approximations to slow 
rotator ellipticals, which are weakly triaxial and have spherical 
isophotes in the inner half-light radius (Simon et al. 2024). 
Cylindrically aligned velocity ellipsoids are typically applied to fast 
rotator massive ellipticals (Cappellari 2008) but have also been 
assumed for slow and nonregular rotators, as in the ATLAS3d 
survey (Cappellari et al. 2013). Here, we compare JAM models 
using both a cylindrically aligned and spherically aligned velocity 
ellipsoid. In the cylindrical case, we assume radial anisotropy 
B-—1-— a2fa5, is constant with c7, — bo?. In the spherical 
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Table 4 
The Self-consistent Model Posterior Medians and 1c Ranges for the Different Kinematic Models/Data Sets 
S/N Velocity Ellipsoid Central Kin. Bin i (deg) 69/0, o/o M. (M5) Y (Mo/Lo,i) 
20 Spherical Kept 731 1.0703 e Q.8*08) x 1010 4.9*02 
20 Spherical Removed 7S 0.83013 e (2.0*09) x 1010 405 
20 Cylindrical Kept 66*13 - 0.927019 (2.7703) x 1079 4.8503 
20 Cylindrical Removed 66*17 I: 0.91*048 (2.4703) x 1010 4.9703 
18 Spherical Kept 72312 0.85010 e (2.0*09) x 1010 44107 
18 Spherical Removed 7213 0.837018 e (1.9488) x 10!° 4.1738 
18 Cylindrical Kept 66*13 : 0.89*045 (2.3*03) x 1010 4204 
18 Cylindrical Removed 68*12 0.897933 (2.2704) x 1010 424 
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Figure 4. Vms = yV? + o? vs. radius compared with the best-fit spherical JAM $ 
model. While the fit was to all data, the points are color-coded by major axis/minor > 
axis sector (Figure 1), and the curves show the expected model runs along the major 
and minor axes. The V and c values are available. ee 96 OS OP AD o D a^ Oo a^ S 2 
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case, radial anisotropy 3 = 1 — 07/0? is constant. The JAM 
models include no tangential anisotropy y = 1 — 0; / o2 = 0, 
since the velocity second-moment predictions are independent of 
"y, nonzero tangential anisotropy is not usually needed to model the 
overall kinematics of real galaxies (Cappellari 2016). 


3.1. JAM Models 


For both the spherically and cylindrically aligned cases, we 
incorporate a self-consistent model, assuming that the total 
mass density follows the surface brightness profile. In 
particular, the projected mass density can be obtained by 
multiplying each Gaussian component in Equation (1) by a 
constant mass-to-light ratio Y». The parameters for this model 
are then Yy, the radial anisotropy parameter (c;/og for 
cylindrical or c/o, for spherical), inclination i, and the black 
hole mass M.. Here inclination i is the angle between the 
galaxy's symmetry axis and the line-of-sight direction and M. 
represents the total mass of the binary, which is modeled by a 
Keplerian point mass potential. We represent the PSF of the 
kinematic observations with a single circular Gaussian of 
dispersion 0727 (as described in Section 2). We do not 
incorporate a model for the contribution of the dark halo's 
potential due to the limited S/N and radial range of the 
kinematic data; tests with a Navarro—Frenk-White halo give no 


Figure 5. Corner plot for the spherical S/N — 20 JAM fit. 


halo constraints and negligible changes to the other fit 
parameters. 
At S/N ~20 the central bin has a large Vins = 


JV? + o? = 532 kms! (Figure 4). Since this has a strong 
effect on the black hole mass fit, we tested for sensitivity and 
systematic bias by optionally dropping the innermost bin. We 
also fit using a lower Voronoi target S/N — 18, for additional 
bins across the central 0/5. In all cases, the V, rise from the 
SMBHB is still apparent. 


3.2. Results 


We perform Markov Chain Monte Carlo fitting of the JAM 
models, with a sin(7) prior on the inclination and wide uniform 
priors on the other parameters. The marginal posterior values 
and 6846 confidence regions are summarized in Table 4 for two 
S/N binnings, both velocity ellipsoid models, and retention or 
exclusion of the central kinematic point. We note that while 
error bars are larger on the spherical model parameters and 
while low S/N binning tends to slightly reduce the best-fit 
black hole masses, all results are consistent at the ~1ø level. 

Some additional systematic sensitivity comes from the 
choice of the spectrum/template polynomials. We explored 
this by refitting the cylindrical S/N — 20 model (with its 
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Figure 6. (a) Correlation between core radius and black hole mass from Dullo et al. (2021), including measurement of depleted core Holm 15A from Mehrgan et al. 
(2019). (b) Correlation with effective velocity dispersion for BCGs from Bogdan et al. (2018), represented by the orange data points and blue fit line. The black line is 
the trend for all core ellipticals from Saglia et al. (2016). Note that 4C+37.11/Holm 15 lie well above the relation for all ellipticals, but are consistent with the BCG 


trend. 


smaller statistical errors), varying the degree of multiplicative 
and additive polynomials from 4 through 7 (16 combinations). 
Since the SMBHB mass M. is our prime measured parameter, 
we focus on its sensitivity. While over half of the values fell 
within 0.5 x 10'?M... of the 4A+6M fit result, values ranged 
from 2.25 — 3.75 x 10'?M... Evidently, significant systematic 
uncertainties remain in the velocity extractions; these can be 
controlled with higher S/N data. 

Since 4C--37.11 appears to be a slow rotator, we take the 
spherical S/N —20 fits as our standard, giving M.— 
(2.8 + 0.8) x 10? Mz, statistical errors only. While the mass 
is large, so are the error bars, encompassing nearly all of the 
other fits and polynomial choices. However, a slightly more 
conservative interpretation might adopt M. 72.4 x 10? M... 
about 0/2 lower, since the central S/N = 20 bin does tend to 
increase the mass. The radial V,,,, run and corner plots for the 
spherical fit are shown in Figures 4 and 5. In the latter, we note 
that i is only weakly constrained within the prior and that M. is 
correlated with the velocity anisotropy, with larger hole masses 
demanded for 3 « 0. Higher S/N is clearly needed to pin down 
B (and its possible radial dependence due to core scouring, see 
below). 

Using the black hole mass, constant stellar mass-to-light ratio, 
and inclination for the spherical model, we estimate the radius of 
the sphere of influence rsor of the binary. Following Mehrgan 
et al. (2019) and Simon et al. (2024), we define rsoi by 
M. & rsop) =M., where M, is the enclosed stellar mass. Using 
the routine mge radial mass of JamPy (Cappellari 2008) to 
compute the enclosed stellar mass within a radius r using the 
MGE of the surface brightness, we find rgo; 1/29 (1.38 kpc), 
quite consistent with the observed Vms flattening. 


4. Discussion and Conclusions 


It is interesting to compare our dynamical M.~ 
2.8 x10'°Mo estimate for 4C+37.11’s binary with values 
found for other massive ellipticals. These M. values correlate 
with host properties, especially those of the host core, giving 
insight into black holes/host coevolution (Kormendy & 
Ho 2013; Saglia et al. 2016; Bogdán et al. 2018). 


The core radius 7, has been found to tightly correlate to the 
central black hole mass among massive ellipticals (Mehrgan 
et al. 2019; Dullo et al. 2021). We compare our 4C4-37.11 
values with dynamically measured M. and r, for other elliptical 
galaxies (Dullo et al. 2021), including Holm 15A, which has 
the largest dynamically measured black hole mass in the local 
Universe, from Mehrgan et al. (2019). We find our value is 
quite consistent with these authors' best-fit trends for cored 
ellipticals (Figure 6(a)). 

Another well-known correlation is the M.—o relation, which 
has been found to depend on the host type. While our 
oe = 332+11kms ' is somewhat suspect since we can only 
integrate to ~0.25R,, we can compare with the trends seen by 
various authors. We first consider the cored ellipticals with 
dynamical mass measurements in Saglia et al. (2016). Our 
estimate exceeds their best-fit correlation’s prediction by 
~7.2 x. However, if we compare with the cored ellipticals 
that are also brightest cluster galaxies (BCGs) from Bogdan 
et al. (2018), we find that 4C+37.11 (and Holm 15A) follow 
their best-fit trend more closely (Figure 6(b)). 

The M.~ 2.8 x 10'°M., dynamical mass for 4C+37.11’s 
SMBHB is one of the largest measured in the local Universe. 
While exceeding the value predicted for a typical elliptical, it is 
in line with expectations for cored ellipticals and, especially, 
those that are also BCGs. This fits well with the picture that the 
host is a fossil cluster, the product of several major mergers. If 
these mergers were largely dry, dissipation-less events, they 
would grow the central black hole mass faster than the stellar 
velocity dispersion (Mehrgan et al. 2019). In addition, the 
back-action from the current binary (and likely from earlier 
binary phases) has “scoured” the core, removing stars capable 
of exerting dynamical friction via three-body gravitational 
slingshots (Begelman et al. 1980). In that, 4C--37.11 resembles 
other extreme cored ellipticals, including Holm 15A. This 
picture can be tested and extended by additional kinematic 
studies, which can measure /X(r), tighten the M. estimate, and 
constrain the orbit anisotropy of the stars in the host core, 
probing the system's merger history. While such measurements 
will be challenging or impossible from the ground, the JWST 
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NIRSpec IFU field is well-matched to the host r, and can 
provide important information on the core kinematics. 
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